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The electronic properties of one-dimensional graphene superlattices strongly depend on the atomic size and 
orientation of the ID external periodic potential. Using a tight-binding approach, we show that the armchair 
and zigzag directions in these superlattices have a different impact on the renormalization of the anisotropic 
velocity of the charge carriers. For symmetric potential barriers, the velocity perpendicular to the barrier is 
modified for the armchair direction while remaining unchanged in the zigzag case. For asymmetric barriers, 
the initial symmetry between the forward and backward momentum with respect to the Dirac cone symmetry 
is broken for the velocity perpendicular (armchair case) or parallel (zigzag case) to the barriers. At last, Dirac 
cone multiplication at the charge neutrality point occurs only for the zigzag geometry. In contrast, band gaps 
appear in the electronic structure of the graphene superlattice with barrier in the armchair direction. 

PACS numbers: 72.80.Vp 73.22.Pr 


I. INTRODUCTION 

Graphene has attracted much attention since 20041 with 
the first experimental characterizations unveiling its amaz¬ 
ing electronic and transport properties!. The unique behav¬ 
ior of massless chiral electrons present in graphene gives rise 
to uncommon effects such as the minimum of conductiv¬ 
ity!, Klein tunneling! of charge carriers and the anomalous 
quantum hall effect!. Such tremendous properties related to 
pseudo-relativistic effects have motivated intense efforts to¬ 
wards the realization of graphene-based electronics! - —. 

Within the variety of graphene-based devices for which a 
great potential is foreseen, the specificities of graphene su¬ 
perlattices have triggered special interest!! - —. In particular, 
one-dimensional (ID) graphene superlattices (SLs), which 
can be obtained by controlled surface patterning, might ex¬ 
hibit electron beam supercollimation!! Furthermore, the 
possibility of creating new Dirac cones at the charge neu¬ 
trality point (CNP) or at low excitation energy is also an 
interesting feature of graphene SLs, both in single and bi¬ 
layer graphene!!. The appearance of new Dirac cones cor¬ 
responds to observable dips in the density of state (DOS), 
which were observed experimentally by STM measurements 
on graphene exposed to a Moire pattern!! - — and in corru¬ 
gated graphene!!. ID SLs can serve as simplified models 
for both corrugated graphene and graphene modulated by an 
external electrical field. 

The directional dependency of the electron propagation in 
ID graphene SLs is related to an unexpected anisotropy of 
the electron velocity. In the direction perpendicular to the 
potential barrier, the velocity of the electron remains con¬ 


stant and is always equal to the velocity of pristine graphene. 
This behavior is unaffected by the width (W), the height ( U ) 
or the modulation period (L) of the barrier. In sharp con¬ 
trast, the velocity parallel to the barrier is strongly impacted 
by these three parameters. In certain cases, the velocity in 
the direction parallel to the barrier can even completely van- 
ish!!*!!. Further experimental work confirmed and extended 
these first observations. In particular, Dubey et al. JH suc¬ 
ceeded to fabricate ID SLs and observed the variation of the 
resistance of this structure as a function of the barrier height. 
Their observations corroborate the appearance of new Dirac 
cones at the charge neutrality points predicted by Barber et 
al (32] and Ho et al (33tl . Including the influence of the bar¬ 
riers edges, Lee et al [34|] showed by ab-initio calculation 
that gaps can appear in ID SLs with barriers in the armchair 
direction, while this is not the case in the zigzag direction. In 
the experimental part, the graphene samples were suspended 
over periodic nanotrenches patterned in the substrate. Since 
interaction with the substrate is believed to introduce chem¬ 
ical doping through the adsorption of oxygen, the periodic 
potential was modeled by epoxy oxygens ad-atoms which 
can lead to impurity specific physics, such as resonant states 
leading to localization effects!!. It is therefore not clear if 
the the apparition of those gaps are due to the oxygens ad¬ 
atoms or to the presence of the SL. 

The first results by Lee et al 0 and the known im¬ 
portance of edge physics in graphene ribbons!!*!! call for 
a more systematic analysis, taking into account the effect of 
edge physics of the barriers in ID SLs. The experimental 
realization of controlled ID SLs is possible, albeit a chal¬ 
lenging task. If the chiral directions of a graphene sheet are 
known, alignment of periodic gates on a patterned surface is 
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possible. As an example, Ponomarenko et al. successfully 
align graphene on hexagonal boron nitride^ . The knowl¬ 
edge of the chirality is more challenging. A first possibility 
is through CVD grown graphene which was reported to ex¬ 
hibit hexagonally shaped grains with aligned edges, mostly, 
in the zigzag direction^. Another solution would be to etch 
a small portion of the graphene via iron catalysts^. As the 
chirality along the etched line is preserved, alignment be¬ 
comes possible in this direction. 

In this article, the differences between the zigzag SLs 
(ZSLs) and armchair SLs (ASLs) are investigated system¬ 
atically by comparing the impact of the three relevant pa¬ 
rameters defining a ID SL: E7, W and L. The following 
conclusions are drawn. By taking into account the different 
orientations of the ID SL (zigzag and armchair), the Dirac 
cone multiplication behavior is complexified. The energy 
position of the additional Dirac cones depends on the value 
of U, as predicted by Park et al. in. In addition, the new 
cones at the charge neutrality point (CNP) described by Ho 
et al. L33] only appear for ZSLs and lead to the apparition 
of three families of cones. In contrast, for ASLs, band gaps 
appear in the density of states. Depending on the number of 
carbon dimers making up the barrier, ASLs can be classified 
into three families. Finally, the velocity modulation quanti¬ 
fied by Barbier et al El is also impacted by the SL barrier 
edge geometry. For symmetric barriers, the velocity perpen¬ 
dicular to the barriers is modified for ASLs and is robust in 
the case of ZSLs. For asymmetric barriers, velocity renor¬ 
malization is observed in both directions for all SLs. These 
asymmetric barriers also break the initial symmetry between 
the forward and backward momentum direction for one of 
the two main directions of the velocity: in the perpendicular 
direction for ASLs and in the parallel direction for ZSLs. 

This systematic analysis is organized into three different 
Sections in this paper, namely the impact of the parameters 
defining the barrier on the DOS and the band structure (Sec- 
tion[III]), on the velocity at the CNP (Section[IV]) and finally, 
on the velocity at higher energies (Section[V]). The numerical 
techniques and the model used to describe the SLs are pre¬ 
sented in Section ITT1 Section PvTl describes the effect of white 
noise disorder on our calculations. 


II. NUMERICAL TECHNIQUES AND MODEL 

The Hamiltonian is expressed in an orthogonal single p z 
orbital basis set. The tight-binding model accounts only for 
first nearest-neighbors interactions described by the hopping 
term 70 =-2.6 eV. The Hamiltonian of the two-dimensional 
periodic system reads as 

H = € z Y J 4 C i+"Y0 J2 c\cj ( 1 ) 

i <i,j> 

where c\ and <7 are creator and annihilator operator on 
atomic site i and where < i,j > denotes the sum runs only 


(a) Zigzag (b) Armchair 




FIG. 1: Representation of graphene SLs with the potential applied 
in (a) the zigzag direction and (b) the armchair direction. The blue 
regions contain atoms with an applied potential a that form poten¬ 
tial barriers of width W\. The red regions form barriers of width 
W 2 and contain atoms with applied potential 62 . L is the periodic¬ 
ity of the pattern. 


on atomic site j being first nearest neighbors of atomic site 
i. The onsite energy term (e z ) accounts for the local electro¬ 
static environment. 

To create a ID periodic superlattice potential on top of the 
2D graphene plane, two electrostatic regions can be defined 
by setting e z to e\ or 62 depending on whether the atomic 
site belongs to the first or the second electrostatic region. No 
other modifications than the onsite terms of the Hamiltonian 
are performed to account for the presence of the electrostatic 
barrier potential. In particular, the hopping terms between 
atomic sites belonging to the two separate regions remain 
unchanged (i.e. equal to 70 ). While cutting graphene into 
nanoribbons will produce physical edges, here the system 
remains an infinite 2D graphene plane but with a ID periodic 
potential imposing a given ID orientation with respect to the 
crystal. 

Figure Q] depicts how ZSLs and ASLs are modeled. Their 
structures are composed of two potential barriers of widths 
W\ and W 2 repeated periodically with period L = W\ + 
W 2 . The barrier heights are given by and 62 respectively. 
Actually, only the potential difference U = e\ — 62 

induces a modification of the electronic properties dis¬ 
cussed in this paper, notwithstanding a rigid shift of the CNP, 
which is given by 

W 1 e 1 + W 2 e 2 ^ 

L • ( j 

It is therefore easier to work with an effective potential 
U by setting e\ = 0 eV and e 2 = U, and to realign a pos¬ 
teriori the CNP to zero energy in order to compare the dif¬ 
ferent systems. Only moderate values of U are considered, 
as larger values are not achievable experimentally and be¬ 
cause the model becomes unreliable for large value of U. 
Finally, in this article, W = W 2 and, unless stated other¬ 
wise, W 2 = L /2 (symmetric barrier). 

The presence of a ID periodic potential, creating a regular 
superlattice, maintains the transport in the ballistic regime, 
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i.e. the propagation is governed by the periodic extended 
Bloch states. Therefore in such perfect SLs, the only rele¬ 
vant quantities governing the ballistic transport are the DOS 
and the carrier velocities. The velocity can be evaluated from 
electronic band structure computed by a direct diagonaliza- 
tion of the Hamiltonian around a given k point: 


1 d E(k) 


Vk h dk 


(3) 


The energy-dependency of the velocity v(E) can then be 
computed by integrating over the whole Brillouin zone. Al¬ 
ternatively, the real-space TB-based Kubo-Greenwood for¬ 
malism, as described in Refs. l3l4lU 47ll, can be used to 
compute carrier velocities from the spreading of wavepack- 
ets in the ballistic regime. Unfortunately for graphene, the 
velocity cannot be accurately computed at the CNP in this 
formalism because a mathematical singularity exists at this 
point (see appendix [A| for more information). The other 
transport quantities such as the semi-classical conductiv¬ 
ity, the mean free path or the mobility, typically accessi¬ 
ble within the Kubo-Greenwood approach, are only well de¬ 
fined into the diffusive regime. This regime is however only 
obtained after multiple scattering events in presence of a 
stochastic disorder potential. The Kubo-Greenwood method 
can also describe quantum localization effects in disordered 
systems. The present paper mainly focuses on perfect sys¬ 
tems, where charge carriers remain in the ballistic regime 
for the whole energy spectrum. The impact of disorder is 
discussed in Section [VT] at the end of the paper. 

Finally, the DOS can be calculated independently of the 
time evolution of wavepackets thanks to the Hay dock recur¬ 
sion method using a Lanczos algorithm with continued frac¬ 
tion calculations^. 


III. DENSITY OF STATE AND ELECTRONIC 
BANDSTRUCTURE 


pertinence there and because such energies are in principle 
experimentally inaccessible to electronic transport measure¬ 
ments. 

In the DOS, each dip associated to a Dirac cone is pre¬ 
ceded by a peak due to secondary Van Hove singularities 
(VHS) in the electronic band structure. The amplitude of 
these VHS is influenced by the potential U. For small val¬ 
ues of U, the dips and associated secondary VHS are not 
visible in the DOS although Dirac cones can be identified in 
the electronic band structure (not shown here). 

For values of W different from the symmetric case (i.e. 
W 7 ^ L/2), an electron-hole asymmetry in the DOS is cre¬ 
ated, and the amplitudes of the secondary VHS change. Nev¬ 
ertheless, the energy position of the peaks is robust to this 
change of W. For a given L, the DOS curves corresponding 
to SLs given by widths W and |L/2 — W\ are antisymmet¬ 
ric to each other around the CNP [see panels (b) and (c) of 

Fig. El 



Energy (eV) Energy (eV) 


A. Generic characteristics for ASLs and ZSLs 

The main impact of a SL on the DOS is the presence of 
the local dips corresponding to the apparition of new Dirac 
cones. The energy position of these dips in the calculated 
DOS presented in Fig. [2] using U = 1.04 eV (red dots for 
ASL and blue line for ZSL) is in agreement with the follow¬ 
ing analytical formula (vertical dashed lines) given by Park 
et al. lfl2h up to E = 1 eV , which confirms that the position 
of the new Dirac cones only depends on L: 

777 7T 

E m = ±hv F -j- (4) 

with v f the Fermi velocity in pristine graphene and m an in¬ 
teger. At higher energies, the agreement between numerics 
and analytics becomes gradually worse, but this is not rele¬ 
vant because the simplified first neighbor TB model loses its 


FIG. 2: DOS of SLs. (a) Difference between the DOS for a ASL 
(dotted red points) and ZSL (solid blue line). Both DOS are given 
for a potential U — 1.04 eV and a period L of 10 nm. The dashed 
vertical lines gives the position of the new Dirac cones as given by 
Eq. (??). The only visible difference is the splitting of the peaks 
for ZSL. Antisymmetric DOS around the CNP for ZSLs with W = 
0.29L (b) and W = 0.71L (c) (L = 10 nm, U = 0.52 eV). 


B. Superlattice barrier edge-dependent characteristics 

For sufficiently large values of U and small values of L, a 
splitting of the secondary VHS occurs for ZSLs. This split¬ 
ting is not observed for ASLs, as depicted in Fig. E^a). In 
addition, for ASLs, band gaps can appear under certain con¬ 
ditions (see Figs. [3][4]and[5]). Their existence depends on the 
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FIG. 3: Band structure close to the K point for two ASLs (U = 
0.845 and U m 2.08, depicted in black and blue, respectively) 
with band gap in the k±( a) and k\\ (b) directions. The CNP was 
not realigned at E — 0 for this figure. As expected, its position is 
given by U/2 (because L = W /2 in Eq. ??). (c) Brillouin zone 
(BZ) associated to the unit cell, the position of the Dirac cone (if it 
exists) is labeled K in reference to its position in the original hexag¬ 
onal BZ of graphene. For ASLs, if multiples cones do not appear, 
the conduction and valance band still present multiple maxima and 
minima in the k\\ direction when L and U are large enough [e.g. 
(b), blue case]. 


number of carbon dimers 3 p + n (with n = 0,1 or 2 and p an 
integer) in the barrier which allows to group the ASLs into 
three families defined by the value of n (see inset of Fig. |4| 
for visual representation of the dimers). At low potential (for 
U > 0), a gap only opens up for the family where n = 0. 
Higher values of potential (U 0) are required to observe 
gaps for the two other families (n = 1,2). The values of 
the gaps are generally very small (few meV) as illustrated in 
Fig. [3 and are inversely proportional to W. The exact value 
of the gap can then be fine-tuned by varying U following the 
bell shape as pictured in Fig. |4[ The DOS corresponding to 
the largest calculated band gap is shown in Fig. [3 


IV. VELOCITY AT THE CHARGE NEUTRALITY POINT 



FIG. 4: Variation of the gap with U for different ASLs. The value 
of W (= L/ 2) is chosen so that all the SLs are from the 3 p family 
(n — 0). The value of L is given in number of dimers as described 
by the inset. 



FIG. 5: Density of state for the ASLs having the largest gap in 
figure 0| (i.e. L containing 6 dimers and U = 2.34eV). The gap is 
clearly visible and has the same value as the one obtained directly 
from de band structure, confirming the existence of the observed 
gaps. 


behavior is discussed in section [V| focusing on the energy 
dependency of the velocity. Large variations of W from the 
symmetric case (W = L/2) reduce the number of cones. In 
other words, if W tends to L or 0 all the new cones disap¬ 
pear. The present section focuses solely on the symmetric 
configuration (W = L/2) where all the new cones appear¬ 
ing in ZSLs are at the CNP. 


Focusing on the central dip in the DOS, i.e. at the CNP, 
the multiplication of cones at this energy observed by Ho et 
al. ll33tl only occurs for ZSLs for which W = L/2. The 
present simulations agree with the creation-by-pairs model 
developed by Barbier et al. 0. Nevertheless, our simula¬ 
tions indicate that these new cones can be classified into two 
categories, each with different properties from the original 
cone (labeled main Dirac cone in the rest of this article). If 
W / L/2, these new cones are shifted in energy away from 
the CNP (positioned at zero energy by convention). This 


A. Symmetric ZSLs 

As predicted by Park et al. El , SLs induce an anisotropic 
velocity renormalization. This picture is depicted and ex¬ 
tended in Fig. [3 The velocities are described as a function 
of U. Because the determining factor for this velocity renor¬ 
malization is actually the product LU , similar curves (not 
shown here) can be obtained by varying L and keeping U 
constant. This comment is valid for the remaining of the 
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FIG. 6: Variation of the velocity at the CNP in the directions perpendicular (_L) (b) and parallel (||) (c) to the potential barrier, for ZSLs. 
Each symbol in (b) and (c) is associated to a particular cone as depicted in panel (a) representing the band structure for a chosen potential, 
U = 2.3 eV, in the direction k\\. L — 10 nm and W — L/2. The main (or original) cone is always at the center and its velocity is depicted 
by the blue circles. The new cones appear two by two, one on each side of the main cone. The last cones to appear are therefore closest to 
the main cone. Corresponding symbols between (a), (b) and (c) panels are used to facilitate reading. 


section. 

In Fig. [6l the main Dirac cone at the CNP (blue circle 
symbol) has the same flavor as the one described by Park et 
al. in. namely that the velocity perpendicular to the bar¬ 
rier, v™ (m for main) in panel (a), is constant and equals the 
velocity in pristine graphene , while the velocity parallel to 
the barrier, v™ in panel (b), varies periodically and goes to 
zero for certain values of U. New cones appear at the CNP 
for symmetric barriers. All the cones generated on the left 
side of the main cone in panel (c) (plus, minus and triangle 
symbols) are part of a second flavor. The third flavor con¬ 
tains the cones created on the right side (diamond, cross and 
star symbols) of the main cone. 

Apart from the two first new cones (plus and star sym¬ 
bols), each set of two new cones appears slightly before the 
minimum of v™ (see, for instance, diamond and triangle 
symbols in Fig. [6^c)). The energy difference between this 
minimum and the energy at which the new cones appear in¬ 
creases with U or L. 

Cones of second and third flavor roughly depict a similar 
behavior for the velocity renormalization with V (or L). In¬ 
deed, v\\ always starts from zero and slowly saturates with 
U ( L). Surprisingly, the velocity for the last cone (star sym¬ 
bols) is higher than the velocity for pristine graphene. The 
reason for it contrasts with the squeezing of the Dirac cone 
due to electron-electron interactions^. A difference in be¬ 
havior between the left and right cones is also visible: the 
velocity of the left cones saturates more quickly than the ve¬ 


locity of the right cones. 

In the perpendicular direction, all characteristics are in¬ 
verted: the velocity goes from the velocity of pristine 
graphene toward zero and decreases more quickly for the 
cones on the right. Barber et al. ll32h found a similar behav¬ 
ior for the additional cones. Nevertheless this separation into 
two classes, with slightly different properties, was missing in 
their analysis. They also found that the new cones appear at 
a minimum of v™, a conclusion which is slightly modified 
here. 

For the velocities calculated at directions in between the 
perpendicular and the parallel one, the velocities v± and v\\ 
always appear for all families as extrema and the velocity 
changes smoothly between them. 


B. Symmetric ASLs 

For ASLs, the impact of the parameters L and U is not 
completely captured by the product LU because of the ex¬ 
istence of three families. When U does not allow to switch 
between those families, the value of L can, and this will thus 
induce a different behavior. 

Figure [7] shows the velocity renormalization for the ASLs 
case (for W = L/2) when changing the value of L and 
figure [8] when the value of V is modified. For the renor¬ 
malization in function of L, all data points are gathered in 
Fig. |7] (a) but also separated into the three families n — 0, 
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L(nm) L(nm) 

FIG. 7: Variation of the velocity with L at the CNP for symmetric 
ASLs for the three families : 3 p, 3p+1 and 3p + 2, with W — L/2 
and U — leV. (a) The datapoints obtained without distinction 
between the families of ASLs depict non-monotonous variations in 
the velocity in the A direction. Separating into the three families 
(b), (c) and (d), the curves make more sense for this direction. In the 
|| direction, the velocity is completely similar to what is observed 
with ZSLs. 


1 or 2 in the panels (b),(c) and (d), according to the previ¬ 
ous discussion in Section IIIIB1 After separation into three 
families, the similarities between variation of U and L are 
recover [compare Fig. 0b),(c),(d) and Fig. 0a),(b),(c) re¬ 
spectively]. The overall behavior of v± is similar as the 
one of ZSLs, provided the values of V or L which create 
a band gap are excluded (inducing an absence of velocity at 
the CNP). Since Dirac cone multiplication does not occur 
for this barrier edge geometry, only the velocity of the main 
Dirac cone is depicted for ASLs. For uy, the behavior is 
similar to ZSLs [compare for instance blue circle symbols in 
Fig. 0c) and Fig. 0 a) (or Figj8])], which confirms the L and 
U -dependency highlighted above. For v±, the behavior is 
different than for ZSLs. Indeed, this velocity now varies, in 
contrast to the constant value observed for ZSLs [Fig. Ob)]. 
As displayed in 0 these variations depend on the family to 
which the ASLs belong to (n = 0,1 or 2), which explains the 
non-monotonous variations in the velocity observed. When 
approaching the gaps, oscillations of the otherwise constant 
value of v± are observed. Perpendicular velocity eventually 
vanishes when reaching the gap region (see FigJS]). 


C. Asymmetric ZSLs 

An asymmetry between the modified and pristine zone 
(W / L/2) has also a large impact on the velocity renormal¬ 
ization. The effect at the CNP is larger for ZSLs. Because of 
the asymmetry, the new Dirac cones are not generated at the 
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FIG. 8: Variation of the velocity with U at the CNP for symmetric 
ASLs for the three families: 3 p, 3p + 1 and 3p + 2. The period L 
used is 21, 19 and 20 dimers respectively. The shaded area high¬ 
lights the position of the gaps. Only the 3 p family has a gap at low 
potential. Near a gap, the oscillations of v± increase up to a point 
where the velocity falls rapidly to zero. 


CNP anymore (and those will thus be discussed in the next 
Section). In this paragraph, only the impact on the main cone 
(blue circle symbols in Fig. I2> is considered. 

The first effect of this asymmetry in ZSLs is to break the 
equivalence between the fey" and fcjj" direction for uy (see 
Fig 02] for sign convention). This difference in the behavior 
of velocities u^ and —u^ is illustrated in Fig. [9j 

The blue circles in panels Fig. 0b) and (c) suggest that, 
for a given valley K, some asymmetry is obtained depend¬ 
ing on the direction of the carriers flow through the barrier 
for U > 1 eV. Nevertheless, a correspondence exists be¬ 
tween the two velocities. The velocity ujj~(—ujj~) obtained 

for a width W is exactly the same as the velocity —u^ (up 
obtained for a width L — W, respectively. 

Looking at the impact of U, both un and v± are affected, 
in opposition to the symmetric case. v± now presents a 
decreasing behavior with U, which gets more pronounced 
when W tends towards L or 0 and completely disappears at 
W = L/2, consistent with previous observations. 

In Fig. 0a), for W = L/ 4, the decrease is not very pro¬ 
nounced in comparison with the (constant) dotted line found 
for W = L/2. As W decreases further, the renormaliza¬ 
tion of v± gets more pronounced [see panel (d)]. The effect 
of U has a larger impact on uy [panels (b) and (c)]. First, 
the minima of uy are not any more at zero. Then, the be¬ 
havior of u^ and u^ is opposite: if the value of the minima 
increases for —u^ [as in panels (c) and (f)], the value of the 
minima in u^ decreases [as in panels (b) and (e)]. In a first 
approximation, if W is not too small, the minima are aligned 
on the red line starting from zero [panels (b) and (c)]. The 
slope of this line is opposite for u^ and —u^ and decreases 
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FIG. 9: Variations of the velocity with U for asymmetric ZSLs, L — 48 dimers 10 nm). The correspondence of the symbols is the same 
as the one used in Fig. [6| (a), (b) and (c): W = L/4 = 12 dimers, (d), (e) and (f): W = L /12 = 4 dimers. When W goes towards 0 (L), 
the asymmetry between the v jj" and vjj" increases. The black arrows in (d), (e), and (f) indicate the direction of the curve when varying W 
towards zero. The reduction is non-linear and increases when W is closest to 0. Similarly, the [/-dependency of the curves is also dilated 
when W decreases (transition between upper and lower panels). At the limit W — 0, the dilation would become infinite, recovering the 
properties of pristine graphene. 


when W goes to L/2, eventually recovering the behavior of 
the symmetric barriers, where both velocity directions are 
degenerated. However, if W is small (a few percent of the 
value of L), the line joining the minima does not cross the 
origin of the plot [i.e. ( 0 , 0 )] anymore, and is different for 
both directions (not shown here). In such case, the larger 
the velocity renormalization is in one direction (take for in¬ 
stance u7), the smaller it is in the opposite direction (ut), as 
already visible by comparing panels (e) and (f). 

Finally, by varying W (see the green dots corresponding 
to the main cone in Fig.QI)]), the impact of asymmetric bar¬ 
riers becomes apparent as well. In this situation, Figs. flOl b) 
and (c) can be obtained from one another by central sym¬ 
metry around the point (v = 0, W = L/2) (as already 
mentioned at the beginning of this paragraph). The other 
curves in this Fig. (blue and red dots in Fig. [10|) are left 
for Section [Vl where the velocities away from the CNP are 
discussed. 


D. Asymmetric ASLs 

For ASLs, the asymmetry induced by the value of W, has 
a weaker influence on the velocity renormalization, in com¬ 
parison with ZSLs. The values of v\\ stay symmetric (i.e. no 
degeneracy lifting between and v7) for W < L/2 and 
the dependency on the U parameter is similar to the one ob¬ 


served for W = L/2. In contrast, the degeneracy is lifted 
for v ±, albeit quite softly. More specifically, the oscillations 
in the vicinity of the band gaps are different, as depicted in 
Fig. HU for the 3p + 1 family (for the position of the band 
gaps for different families, see Fig. [ 8 ]). 


V. ENERGY DEPENDENCY OF THE VELOCITY 

As mentioned previously, the periodic potential does not 
only modify the velocity at the CNP, but also impacts the 
velocity at higher energies. Those energies are the topic of 
this Section. On the one hand (Section [V Alh for energies 
close to the CNP (< O.leV), only the asymmetric ZSLs 
(W 7 ^ L/2) are discussed, for which the new cones shift 
away from the CNP. ASLs do not induce additional cones, 
and are thus excluded from this discussion. On the other 
hand (Section IV B I) , at higher energies (^ O.leV), the dif¬ 
ferences between ASLs and ZSLs disappear. The focus is 
thus shifted towards the velocities at intermediate incident 
angles. 


A. Asymmetric ZSLs 

Going back to Fig. [3 showing the impact of the asym¬ 
metry on the additional cones (plus, minus, cross, and star 
symbols) as a function of U, the inequivalence between fcjj" 

and fcjj" directions is clearly apparent, as it was already ob- 
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FIG. 10: Variations of the velocity with the barrier width when several cones are present for ZSLs. L = 120 dimers 26 nm). A large 
value of L is used to have more points in the curve. Indeed, an integer number of dimers for W and L is required. Therefore, the larger 
the value of L, the larger the number of points that can be computed. U = 0.52 eV. (a) The velocity in the _L direction is constant for the 
main cone and varies from the velocity of pristine graphene to a minimum situated at W = L/2 for the other cones, (b) and (c) represent 
the values of the velocity in the two inequivalent branches of the cones for the parallel direction. Additional cones are only found in the 
central zone with a L between ~16% and ~84%. The variations of velocity of two new cones are opposite. Contrarily to the main cone, an 
inversion of the velocity is observed for the new cones. Results obtained for smaller value of L are similar. 



FIG. 11: Variations of the velocity with U for asymmetric ASLs. 
L — 40 dimers (~ 10 nm) and W = L /10 = 4 (3p+l family). 
No anisotropy exists for vjj". In opposition to ZSLs, v± exhibits an 
anisotropy between the + and — directions. 


served for the main cone (blue circles). 

Varying now the parameter W (Fig. [Tot . two zones can 
be distinguished (separated by dashed vertical purple lines). 
In the first zone (outer region, where W is close to 0 or L ), 
the value of \L — W\ is too small to create multiple cones 
(only the main cone exists). A simple renormalization of 
the velocity parallel to the barrier (up, similar to the one 
by changing the U parameter, is observed. In the second 
zone (central region), the velocities corresponding to the ad¬ 
ditional cones show a more exotic behavior. In particular, 
close to the boundary between the two zones, the velocities 
of the additional cones have the same sign in and up 
This is somehow unusual because it means that the slope of 
the two bands forming the cone in the || direction have the 
same sign. This curious behavior of the additional cones can 
be visualized by looking at the evolution of the band struc¬ 


ture as a function of W in Fig. [121 For W = L/b [panel 
(a)], the cones are tilted in such a way that the velocities have 
the same sign (upper insets zoom in on this peculiar behav¬ 
ior). For higher values of W/L (here at 26.6%), one of the 
branches becomes parallel to the fey axis and thus the veloc¬ 
ity drops to zero [see panel (b)]. After this transition value, 
a situation occurs where the velocities in the cone have op¬ 
posite signs [see panel (c)]. The absolute values of these 
velocities are not yet equal because the cone is still slightly 
rotated. By increasing further the value of W, the velocities 
become closer and closer to each other in absolute value, to 
finally recover the symmetric barriers case discussed in pre¬ 
vious Sections. 


B. Angle dependency 

Up to now, only the 0° (A) and the 90° (||) cases have 
been discussed, being the most straightforward. We found 
that the velocity changes smoothly between these two ex¬ 
treme cases. The present Section extends the analysis for 
these intermediate angles at energies away from the CNP, 
leading to richer physics. For this analysis the velocities are 
computed within the Kubo-Greenwood method in the ballis¬ 
tic regime. ZSLs and ASLs behave very similarly at higher 
energies. Therefore, only ZSLs are depicted. 

In Fig. |T3l both the angle and the energy dependency are 
clearly visible. At the CNP, the evaluation of velocity is hin¬ 
dered by numerical divergences (shaded region in Fig. fl3l) . 
difficult to resolve using the diagonalization trick, because of 
the angle dependency. For energies very close to the CNP, 
the low-angle velocity (close to 0° = v±) is larger than the 
high-angle velocity (close to 90° = v\\). For higher ener¬ 
gies, i.e. away from the CNP, the behavior is inverted. The 
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FIG. 12: Rotation of the Dirac cones with the modification of W for ZSLs. The arrows outside the circle above the graph indicate the 
Active rotation direction when transitioning the values of W from (a) to (c). The arrow inside the cone has no physical meaning, but helps 
visualize the cone rotation.. v + and v~ are defined as the velocity on the two branches of the cones. The corresponding branches are 
color-coded in the mainframes below. v~ is taken on the branch showing an overall decreasing behavior and u + an overall increasing 
behavior, (a) Both branches of the cone have the same slope: v~ and v + have the same sign, (b) One of the bands is flat: v~ or v + is zero, 
(c) The two branches have opposite slopes but the cone is still slightly turned: v~ and v + have opposite signs but different values. 



FIG. 13: Energy dependency of the velocity at different angles, 
for ZSLs. The first curve on the top is for pristine graphene. The 
direction perpendicular to barrier is taken as 0° reference. From 
top to bottom the angles are 90° (uy), 60°, 40° and 0° (u_l) . The 
vertical dashed lines represent the position of the new Dirac points 
generated by the SLs as given by Eq. ??. The potential applied is 
U = 0.52eV with a period of L = lOnm (W = L/ 2). The area 
around the CNP is shaded because the results are numerically too 
unstable there. 


transition between these two regimes occurs around 0.12 eV 
in Fig.[l3l which corresponds to the position of the first peak 
(VHS) in the DOS. 

For even higher energies, periodic oscillations in the ve¬ 
locity appear. The amplitude of these oscillations is maxi¬ 
mum for smallest angles (u_l) and become barely noticeable 
for largest angles. Every minimum of the velocity corre¬ 
sponds to the position of a Dirac cone in the electronic band 
structure, which corresponds also to a minimum in the DOS. 


Even if the minima are not apparent in the DOS (see Fig.0, 
they are clearly visible for the velocity up to approximately 
1 eV in Fig. [13] for low angles. 

As already observed throughout this paper, a minimum in 
the DOS caused by a Dirac point does not necessarily imply 
a maximum in the velocity. This is further confirmed with 
this energy dependent renormalized velocity curve. Chang¬ 
ing from 90° to 0°, a maximum in the velocity can become 
a minimum at the position of the additional Dirac cones. In 
other words, the Dirac cones away from the CNP induce a 
local maximum of the velocity in the direction parallel to the 
barrier and a minimum in the direction perpendicular to the 
barrier. Finally, the energy-dependent oscillations in the ve¬ 
locity are more pronounced with increasing values of U (not 
shown here). 


VI. EFFECT OF DISORDER 

All the potential barriers considered so far were ideal and 
displayed a perfect periodicity, keeping therefore the charge 
carriers in the ballistic regime. The absence of random 
disorder precludes quantum interference phenomena, such 
as weak and strong localization effects. Therefore, to ob¬ 
serve the transport signatures from this paper, experiments 
should aim at minimizing any form of extrinsic disorder 
(ad-atoms, vacancies, trapped screened charged impurities) 
that may lead to strong scattering. In addition, the barriers 
themselves should be free of disorder and atomically 
perfect. This situation is obviously rather far from real 
experimental conditions. Existing literature gives guidance 
on which features should remain robust and which features 
might disappear if such disorder becomes too strong. More 
specifically, for uncorrelated white-noise disorder on the 
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L=10nm 

W=L/2 

U=1.04eV 


pristine 

V=0.0eV 

V=0.5eV 

V=0.8eV 



FIG. 14: Density of states and velocities for ASLs and ZSLs (L =10 
nm, W-Lj 2, 17=1.04 eV) without (V=0 eV) and with Anderson 
disorder (V^=0.5; 0.8 eV). 


barriers, a perpendicular incident angle keeps a robust 
transmission^, suggesting associated transport features to 
be reliable, while transmission is strongly reduced when 
increasing the incident angle. The system should therefore 
be kept as clean as possible to observe parallel conductivity 
features (the perpendicular case being less sensitive). To 
counter the detrimental effect of white-noise disorder on the 
transmission, a long-range correlation between the potential 
barriers may be used2i2i. 

The systematic study of the effect of disorder on elec¬ 
tronic transport in graphene ID SLs is out of the scope 
of this article. However, selected prospective simulations 
of disordered graphene ID SLs are presented in Fig. [141 
The upper panels compare the effect of Anderson (white- 
noise) disorder on the DOS of both ASLs (left) and ZSLs 
(right). Anderson disorder is introduced by randomly vary¬ 
ing the onsite potentials of all atomic sites with a value of 
Ss Pz G [—V/2, +V/2]. The existence and the position of the 
new Dirac points remain noticeably robust, up to V = 0.8 
eV. Equation (??) can still be used to estimate their position. 


However, disorder induces a splitting in the secondary VHS 
for the ASLs, while this splitting already exists in absence of 
disorder for ZSLs. Smoothed curves and VHS peaks may be 
expected by averaging over a manifold of disorder configura¬ 
tions. This makes it very difficult to experimentally differen¬ 
tiate between ASLs and ZSLs, based on the DOS.In Fig. 14, 
both parallel (central panels) and perpendicular (lower pan¬ 
els) velocities are plotted for ASLs (left panels) and ZSLs 
(right panels).The black curves give the velocity for pristine 
graphene, in absence of SLs and Anderson disorder. Apply¬ 
ing the Anderson disorder (yellow and blue curves) on top 
of the clean SL (red curve) does not change much the be¬ 
havior of both velocity directions above ±0.3 eV, which is 
as expected due to the selected V values. Nevertheless, the 
qualitative features (maxima and minima) remain globally 
robust even for low energies. At zero energy, a drop by a 
factor two for V = 0.8 eV is observed for both velocities 
for the ASL case. For the ZSLs case, the changes in ve¬ 
locity at low energy are smaller, with the perpendicular fea¬ 
tures even more robust, in agreement with above statements 
based on literature. Nevertheless, for both SL orientations, 
a large drop (increase) of the perpendicular (parallel) veloc¬ 
ity appears around 0.2 eV, respectively. Further calculations 
would be required to complete this picture, but these prelim¬ 
inary velocity and DOS calculations indicate that ZSLs are 
more robust to the detrimental effect of Anderson disorder 
than ASLs. 


VII. CONCLUSION 


The impact of the ID SL orientation with respect to the 
graphene crystal was examined on grahene ID SLs physics 
were studied. Important differences were highlighted, like 
the presence of new cones at the CNP for zigzag SLs and the 
opening of gaps for armchair SLs. A specific effect induced 
by the SL alignment, absent in the literature, was found in 
the velocity renormalization for the direction perpendicular 
to the ID potential. This renormalization occurs in ASLs, in 
general, and in ZSLs when the barriers are asymmetric. On 
the other hand, the velocity in the direction parallel to the ID 
potential behaves similarly to what is predicted in the litera¬ 
ture, although the position of the minima can be modulated 
by the type of SLs or the parameters used. The asymme¬ 
try of the SLs was shown to have a strong impact on the 
velocity. In particular, it can break the initial symmetry be¬ 
tween the forward and backward momentum direction with 
respect to the Dirac cone symmetry for the velocity in the 
perpendicular direction for ASLs and parallel direction for 
ZSLs. This breaking of the symmetry can be interpreted by 
a rotation and deformation of the Dirac cone(s), leading to 
strong modifications in the associated velocities. By study¬ 
ing the angle dependency of the velocity through the barrier, 
the smooth transition between the parallel and perpendicu¬ 
lar direction is understood. The calculated gaps in ASLs 
are very small. More advanced theoretical frameworks, such 
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as ab initio simulations, are required to correctly assess their 
amplitude. Further studies may focus on systems with mixed 
chiral orientation different than pure armchair or zigzag ori¬ 
entation reproducing certain experimental conditions. Nev¬ 
ertheless, much improved control of edge geometry (using 
a bottom-up approach for chemical synthesis^ - —) make the 
present systems promising for electron collimation experi¬ 
ments. 
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Appendix A: Numerical instability at the charge neutrality 
point 

Numerical instabilities make the accurate calculation of 
the velocity at the CNP difficult. In the following we show 
why and how numerical instabilities appears at the CNR 


This feature is illustrated within the Kubo-Greenwood ap¬ 
proach which we further describe here. 

The Kubo-Greenwood technique was originally devel¬ 
oped to compute the conductivity from the quantum mechan¬ 
ics point of view^2r— and gives information on both the 
quantum and the semi-classical transport. The energy de¬ 
pendent carrier velocities in the ballistic regime can also be 
investigated using this formalism. The formalism is based 
on the propagation of a wavepacket throughout the material, 
described by the diffusion coefficient D defined as: 

m = ^Ai? 2 (i) = X\t) + §- t AY\t) (Al) 

where A R 2 (t) is the mean quadratic spreading of the 
wavepacket. A X 2 (t) and A Y 2 (t) are the mean quadratic 
spreadings in the x and y direction, respectively. The behav¬ 
ior of D with time indicates the transport regime in which 
the wavepacket resides^. If the diffusion coefficient in¬ 
creases linearly with time, the electrons are not scattered 
and move freely in the material (ballistic regime). When this 
coefficient saturates to a certain value (D max ), the electrons 
have experienced sufficient scattering to reach the diffusive 
regime. A further increase or decrease of the coefficient in¬ 
dicates the onset of (anti-)localization^^, rooting in quan¬ 
tum localization corrections. The diffusion coefficient D de¬ 
pends on the mean quadratic spreading of the wavepacket as 

described in Eq. (??), while the spreading in a given direc¬ 
tion (say, x) is given by: 


AX\t) = (\X(t) - X(0)\ : 


Tr 


(xjt) - 1(0 ))* 6(E ~ H) (x(t) - 1 ( 0 )) 


Tr 


S(E - H) 


(A2) 


where the operator X(t) is the position operator in the 
Heisenberg picture. The numerator, and the denominator 
(corresponding to the DOS), are computed separately, with 
the same Lanczos algorithm using continued fractions. To 
reduce the computational cost, the trace is replaced by an av¬ 
erage of random phase states obtained by adding a random 
phase factor to the wavefunction at each orbital of the sys¬ 
tem. Averaging over about 10 random phase states is usually 
enough to reach a satisfactory convergence (< 1% of varia¬ 
tions in the quantities of interest). 

The velocity in the ballistic regime can be extracted from 
the mean quadratic expansion as 


A x\t)=vir^v x = ^Tp± (A3) 


In pristine graphene, at the CNP, both the numerator and 
the denominator of A X 2 {t) [Eq. (??)] tend to zero. The de¬ 
nominator being the graphene DOS, at the Dirac point this 
term obviously tends to zero. From simple physical consid¬ 
erations using Eq. ??, one can show that the numerator of 
A X 2 {t) also tends to zero. Indeed, a finite value of veloc¬ 
ity implies a finite value of A X 2 (t). Since the denominator 
of A X 2 (t) tends to zero in Eq. (??), a finite spreading can 
only exist if the numerator tends to zero too. Mathemati¬ 
cally, using the concept of limits, this is no problem. From a 
numerical point of view, using floating point arithmetics, the 
division of two very small numbers generates large errors, 
explaining the aforementioned instability. A similar prob¬ 
lem can arise when calculating the energy dependent veloc¬ 
ity with a direct diagonalization approach. However in that 
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The Kubo-Greenwood approach was nevertheless helpful 
in this article to compute the energy dependent velocity at 
case, the accuracy is much more controlled and the results an y intermediate angles as shown in Fig [13] 
can be improved by using denser k-point meshes in the Bril- 
louin zone. 
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